/*
 * File:   FriedmannEvolution.cpp
 * Author: gian
 *
 * Created on December 10, 2012, 4:39 PM
 */

#include <math.h>
#include "FriedmannEvolution.hpp"

FriedmannEvolution::FriedmannEvolution(double OmegaRadiation,
        double                                OmegaMatter,
        double                                a0,
        double                                H0) {
  mOmegaRadiation = OmegaRadiation;
  mOmegaMatter    = OmegaMatter;
  mOmegaVacuum    = 1.0 - OmegaMatter - OmegaRadiation;
  ma0             = ma_actual = a0;
  mH0             = H0;
  mEta_actual     = 0.0;
}



void FriedmannEvolution::load(FILE* stream)
{
    fread((void*)&ma0,sizeof(double),1,stream);
    fread((void*)&mH0,sizeof(double),1,stream);
    fread((void*)&mOmegaRadiation,sizeof(double),1,stream);
    fread((void*)&mOmegaMatter,sizeof(double),1,stream);
    fread((void*)&mOmegaVacuum,sizeof(double),1,stream);
    fread((void*)&ma_actual,sizeof(double),1,stream);
    fread((void*)&ma0,sizeof(double),1,stream);
    fread((void*)&ma_previous,sizeof(double),1,stream);
    fread((void*)&mEta_actual,sizeof(double),1,stream);
    fread((void*)&mEta_previous,sizeof(double),1,stream);
}

void FriedmannEvolution::save(FILE* stream)
{
    fwrite((void*)&ma0,sizeof(double),1,stream);
    fwrite((void*)&mH0,sizeof(double),1,stream);
    fwrite((void*)&mOmegaRadiation,sizeof(double),1,stream);
    fwrite((void*)&mOmegaMatter,sizeof(double),1,stream);
    fwrite((void*)&mOmegaVacuum,sizeof(double),1,stream);
    fwrite((void*)&ma_actual,sizeof(double),1,stream);
    fwrite((void*)&ma0,sizeof(double),1,stream);
    fwrite((void*)&ma_previous,sizeof(double),1,stream);
    fwrite((void*)&mEta_actual,sizeof(double),1,stream);
    fwrite((void*)&mEta_previous,sizeof(double),1,stream);    
}

void FriedmannEvolution::set_eta(double eta) {
  mEta_actual = eta;
}

FriedmannEvolution::FriedmannEvolution(const FriedmannEvolution & orig) {}

FriedmannEvolution::~FriedmannEvolution() {}

double FriedmannEvolution::friedmann_f(double a) {
  double r = ma0/a;

  return mH0 * a * sqrt(mOmegaRadiation*r*r*r*r + mOmegaMatter*r*r*r + mOmegaVacuum);
}

void FriedmannEvolution::step(double d_eta) {
  double k1 = d_eta * friedmann_f(ma_actual);
  double k2 = d_eta * friedmann_f(ma_actual + 0.5 * k1);
  double k3 = d_eta * friedmann_f(ma_actual + 0.5 * k2);
  double k4 = d_eta * friedmann_f(ma_actual + k3);

  ma_previous   = ma_actual;
  ma_actual     += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
  mEta_previous = mEta_actual;
  mEta_actual   += d_eta;
}

double FriedmannEvolution::get_OmegaRadiation() {
  double r = ma0/ma_actual;  
  return mOmegaRadiation*r*r*r*r;
}

double FriedmannEvolution::get_OmegaMatter() {
  double r = ma0/ma_actual;

  return mOmegaMatter*r*r*r;
}

double FriedmannEvolution::get_OmegaVacuum() {
  return mOmegaVacuum;
}

double FriedmannEvolution::get_scale_parameter() {
  return ma_actual;
}

double FriedmannEvolution::get_H() {
  return friedmann_f(ma_actual) / ma_actual;
}

double FriedmannEvolution::get_eta() {
  return mEta_actual;
}


//~ Formatted by Jindent --- http://www.jindent.com
